
#include <stdlib.h>
#include <stdio.h>
#include <string.h>
#include <assert.h>
#include <math.h>
#include <sys/time.h>

#include <cutil.h>
#include "mlsKernel.cu"

double get_clock();
void runMLS(int argc, char** argv);

int main(int argc, char** argv) {
   
   runMLS(argc, argv);
   
   //~ CUT_EXIT(argc, argv);
}

void setControlPoints(float4 *P, float4 *Q, const unsigned int N) {
   
   //~ P[0].x =-1.0; P[0].y = 0.0; P[0].z =-1.0;
   //~ P[1].x =-1.0; P[1].y = 0.0; P[1].z = 1.0;
   //~ P[2].x = 1.0; P[2].y = 0.0; P[2].z =-1.0;
   //~ P[3].x = 1.0; P[3].y = 0.0; P[3].z = 1.0;
   //~ P[4].x = 0.0; P[4].y = 0.0; P[4].z = 0.0;
   
   //~ Q[0].x =-1.0; Q[0].y = 0.0; Q[0].z =-1.0;
   //~ Q[1].x =-1.0; Q[1].y = 0.0; Q[1].z = 1.0;
   //~ Q[2].x = 1.0; Q[2].y = 0.0; Q[2].z =-1.0;
   //~ Q[3].x = 1.0; Q[3].y = 0.0; Q[3].z = 1.0;
   //~ Q[4].x = 0.0; Q[4].y = 1.0; Q[4].z = 0.0;
   
   srand(time(NULL));
   for (unsigned i=0; i<N; ++i) {
      P[i].x = rand()/(float)RAND_MAX;
      P[i].y = rand()/(float)RAND_MAX;
      P[i].z = rand()/(float)RAND_MAX;
      Q[i].x = rand()/(float)RAND_MAX;
      Q[i].y = rand()/(float)RAND_MAX;
      Q[i].z = rand()/(float)RAND_MAX;
   }
}

void genMesh(float4 *mesh_pos, unsigned int mesh_width, unsigned int mesh_height) {
   
   unsigned int i=0;
   float deltax = 1.0/((float)mesh_width-1.0);
   float deltay = 1.0/((float)mesh_height-1.0);
   for (unsigned h=0; h<mesh_height; ++h) {
      for (unsigned w=0; w<mesh_width; ++w, ++i) {
         const float x = 2.0*w*deltax-1.0;
         const float y = -1.0;
         const float z = 2.0*h*deltay-1.0;
         mesh_pos[i] = make_float4(x,y,z,0.0);
      }
   }
   
   //~ printf("xxxxxx initial values after \n");
   //~ for (unsigned i=0; i<mesh_width*mesh_height; i++) {
      //~ const float x = mesh_pos[i].x;
      //~ const float y = mesh_pos[i].y;
      //~ const float z = mesh_pos[i].z;
      //~ printf("%d ", i);
      //~ printf(" %f, %f, %f \n", x,y,z);
   //~ }
}

void runMLS(int argc, char** argv) {
   
   assert(argc == 4);
   
   CUT_DEVICE_INIT();
   
   unsigned int mesh_width, mesh_height, NCTRLPOINTS;
   mesh_width = atoi(argv[1]);
   mesh_height = atoi(argv[2]);
   NCTRLPOINTS = atoi(argv[3]);
   
   //~ printf(" mesh_width %d  \n", mesh_width);
   //~ printf(" mesh_height %d  \n", mesh_height);
   //~ printf(" NCTRLPOINTS %d  \n", NCTRLPOINTS);
   
   // allocate host memory for control points
   //~ unsigned int NCTRLPOINTS = 5;
   unsigned int size_P = NCTRLPOINTS;
   unsigned int mem_size_P = sizeof(float4)*size_P;
   float4* h_P = (float4*)malloc(mem_size_P);
   unsigned int size_Q = NCTRLPOINTS;
   unsigned int mem_size_Q = sizeof(float4)*size_Q;
   float4* h_Q = (float4*)malloc(mem_size_Q);
   // initialize control points (in host memory)
   setControlPoints(h_P, h_Q, NCTRLPOINTS);
   
   // allocate device memory for control points
   float4 *d_P, *d_Q;
   CUDA_SAFE_CALL(cudaMalloc((void**) &d_P, mem_size_P));
   CUDA_SAFE_CALL(cudaMalloc((void**) &d_Q, mem_size_Q));
   
   // copy host memory to device
   CUDA_SAFE_CALL(cudaMemcpy(d_P, h_P, mem_size_P, cudaMemcpyHostToDevice));
   CUDA_SAFE_CALL(cudaMemcpy(d_Q, h_Q, mem_size_Q, cudaMemcpyHostToDevice));
   
   // Bind the array to the texture
   cudaChannelFormatDesc channelDesc4 = cudaCreateChannelDesc<float4>();
   //~ cudaChannelFormatDesc channelDesc = cudaCreateChannelDesc(32, 32, 32, 32, cudaChannelFormatKindFloat);
   CUDA_SAFE_CALL(cudaBindTexture(0, texP, d_P, channelDesc4));
   CUDA_SAFE_CALL(cudaBindTexture(0, texQ, d_Q, channelDesc4));
   
   // allocate host memory for vertex positions
   unsigned int nvertices = mesh_width*mesh_height;
   unsigned int mem_nvertices = sizeof(float4)*nvertices;
   float4* h_pos = (float4*)malloc(mem_nvertices);
   // generate vertex positions (in host memory)
   genMesh(h_pos, mesh_width, mesh_height);
   
   // allocate device memory for vertex positions
   float4* d_pos;
   CUDA_SAFE_CALL(cudaMalloc((void**) &d_pos, mem_nvertices));
   
   // copy host memory to device
   //~ CUDA_SAFE_CALL(cudaMemcpy(d_pos, h_pos, mem_nvertices, cudaMemcpyHostToDevice));
   
   cudaArray* cu_array;
   CUDA_SAFE_CALL(cudaMallocArray(&cu_array, &channelDesc4, mesh_width, mesh_height));
   CUDA_SAFE_CALL(cudaMemcpyToArray(cu_array, 0, 0, h_pos, mem_nvertices, cudaMemcpyHostToDevice));
   CUDA_SAFE_CALL(cudaBindTextureToArray(texPositions, cu_array, channelDesc4));
   
   double t_start = get_clock();
   
   // execute the MLS kernel
   dim3 block(8, 8, 1);
   dim3 grid(mesh_width/block.x, mesh_height/block.y, 1);
   MLSkernel<<<grid, block>>>(d_pos, mesh_width, mesh_height, NCTRLPOINTS);
   
   // check if kernel execution generated and error
   CUT_CHECK_ERROR("Kernel execution failed");
   
   // copy result from device to host
   CUDA_SAFE_CALL(cudaMemcpy(h_pos, d_pos, mem_nvertices, cudaMemcpyDeviceToHost));
   
   // stop and destroy timer
   //~ printf("initial values after \n");
   //~ for (unsigned i=0; i<nvertices; i++) {
      //~ const float x = h_pos[i].x;
      //~ const float y = h_pos[i].y;
      //~ const float z = h_pos[i].z;
      //~ printf("%d ", i);
      //~ printf(" %f, %f, %f \n", x,y,z);
   //~ }
   double mlsTime = get_clock() - t_start;
   //~ printf(" mls CUDA Time %f  sec. \n", mlsTime);
   printf("%f \n", mlsTime);
   
   //~ printf("aaaaaaaaaaaaaaaaaa \n");
   
   CUDA_SAFE_CALL(cudaUnbindTexture(texP));
   CUDA_SAFE_CALL(cudaUnbindTexture(texQ));
   CUDA_SAFE_CALL(cudaUnbindTexture(texPositions));
   
   // clean up memory
   free(h_P); 
   free(h_Q); 
   free(h_pos); 
   CUDA_SAFE_CALL(cudaFree(d_P)); 
   CUDA_SAFE_CALL(cudaFree(d_Q)); 
   CUDA_SAFE_CALL(cudaFree(d_pos)); 
}

double get_clock() {
   
   struct timeval t;
   struct timezone tz;
   gettimeofday(&t,&tz);
   return (double) t.tv_sec + (double) t.tv_usec*1E-6;
}
